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Spectral  element  methods  are  high  order  methods  that  combine  the  flexibility  of  finite 
element  methods  with  the  "infinite  order  accuracy"  of  spectral  methods.  The  domain  of  computation 
is  decomposed  into  some  subdomains  -  the  elements  -  (generally  these  are  deformed 
parallelotopes),  and  the  exact  solution  is  approximated  by  a  piecewise  polynomial  of  high  degree. 
The  spectral  element  method  differs  from  other  spectral  methods  using  domain  decomposition 
techniques  (patching  methods)  by  the  way  the  matching  conditions  are  handled.  These  are,  like  in 
the  finite  element  method,  implicitly  taken  into  account  by  the  variationnal  statement  of  the 
discrete  problem.  This  allows  for  more  flexibility  with  no  loss  of  the  spectral  accuracy  (  see,  e.g., 
[P],  [M.P.]  and  [F]).  When  the  algebraic  equations  resulting  from  this  kind  of  discretization  are 
obtained,  the  problem  that  remains  is  to  solve,  in  an  efficient  way,  the  algebraic  system. 

The  interest  of  domain  decomposition  t^hnique  is  to  fraction  the  computational  task  so  as  to 
yield  smaller  problems  and  to  use  parallel  computers  for  instance.  If  the  value  of  the  approximate 
solution  were  known  on  the  various  interfaces,  the  problem  would  be  very  simple  since  it  would 
results  into  the  resolution  of  as  many  disconna^tol  problems  as  the  number  of  elements.  The  main 
difficulty  is  that  this  value  is  not  known;  hence  a  technique  known  as  the  iteration  per  subdomains 
has  been  proposed  in  the  literature  [F.Q.Z.]  [Q.Sl]  to  discover  this  value  iteratively.  Another 
approach  is  to  try  to  invert  the  whole  problem  by  not  working  iteratively  on  each  subdomain.  This 
method,  [M.P.]  (F.R.D.P.]  [R],  consists  in  reducing  iteratively  the  residue  at  the  same  time  on 
every  subdomain.  The  global  method  used  can  be  based  on  a  conjugate  gradient  algorithm  or  another 
iterative  procedure, 

In  the  first  part  of  this  paper,  E.  M.  Ronquist  and  A.  T.  Patera  have  presented  some 
results  concerning  a  new  multigrid  method  for  the  resolution  of  the  algebraic  system  resulting 
from  the  approximation  of  a  second  order  P.O.E.  by  spectral  or  spectral  element  method  in  the 
one-dimensional  domains.  The  very  simple  idea  of  using  the  Jacobi  preconditioner  as  a  smoother 
for  the  iterative  multigrid  algorithm  appears  to  be  a  very  good  one.  Indeed,  the  numerical 
properties  of  this  smoother  seems  to  surpass  all  expectations;  the  reduction  rate  of  each  V-cycle 
appears  numerically  to  be  independent  of  the  discretization  parameters. 
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When  iterative  techniques  are  used,  it  is  important  to  understand  why  these  methods 
converge  in  order  to  foresee  the  generalization  and  the  ability  of  the  methods  to  be  adapted  to  more 
than  test  problems.  Here  we  propose  an  analysis  of  this  phenomenon  and  provide  the  justification 
of  these  very  good  properties.  Many  general  convergence  proofs  exist  in  the  literature  for  the 
numerical  analysis  of  the  multigrid  technique;  among  them  let  us  cite  [M.M.],  [B.D.].  We  use  here 
the  abstract  framework  developed  by  R.  E.  Bank  and  C.  C.  Douglas  [B.D.]  that  fits  exactly  the 
numerical  conclusion  of  [R.P.]  concerning  the  optimal  choice  of  the  smoothing  operations. 

To  our  knowledge,  the  numerical  analysis  of  the  convergence  for  the  multigrid  algorithms 
used  in  spectral  type  techniques  is  somehow  empty.  The  main  reason  is  certainly  that  the  nice 
analysis  that  can  be  done  requires  a  variational  framework  and  the  awareness  that  the  spectral 
methods  are.  exactly  or  very  close  to,  variational  approach  is  not  so  old.  The  other  reason, 
perhaps,  is  that  the  previous  multigrid  techniques  applied  to  spectral  type  methods  [Z.W.H.  1 ,2] 
used  a  finite  difference  preconditioner  as  a  smoother  and  a  Chebyshev  framework.  The 
convergence,  in  this  case,  is  not  so  brillant  as  in  the  present  approach  and  a  priori  more  related  to 
the  good  properties  of  preconditioner  of  the  finite  difference  operator.  Besides  the  variational 
formulation  involves  a  nonsymmetric  form  that  makes  the  analysis  much  more  involved. 

The  paper  is  organized  as  follows;  in  swjtion  II,  we  recall  the  theory  of  [B.D.]  in  a  form 
adapted  to  our  analysis.  In  section  1 1 1 ,  we  first  explain  on  the  test  example  of  the  Galerkin  spectral 
approximation  of  the  homogeneous  Poisson  problem  the  fundamental  reasons  of  optimal  properties 
of  this  multigrid  method.  The  tools  are  based  on  Jackson  inequality  and  some  refined  version  of  the 
approximation  properties  of  the  L^-projection  operator.  In  section  IV,  we  generalize  the  analysis 
to  the  case  of  the  spectral  element  approximation.  We  compare  in  each  section  the  results  obtained 
by  the  theory  with  the  numerical  results  presented  in  [R.P.].  The  last  section  V  deals  with  the  one 
domain  multigrid  technique  when  applied  to  a  non  constant  second  order  problem. 

The  generalization  of  these  results  to  multidimensional  problems  will  be  presented  in  a 


future  paper. 
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ti^P_os1t1on  of  the  problem  and  abstract  theorem. 

11.1  Generalities  on  variational  muUlgrId  techniques. 

In  this  subsection,  we  first  recall  the  theory  developed  by  R.  Bank  and  C.  Douglas  to  analyze 
the  convergence  rate  of  the  muUigrid  algorithm  for  solving  the  linear  algebraic  system  that  arises 
from  the  numerical  approximation  of  elliptic  partial  differential  equations.  We  present  it  in  a 
version  that  we  shall  use  afterwards.  First  of  all,  let  36  be  a  Hilbert  space,  a  be  a  continuous, 
elliptic,  symmetric,  bilinear  form  and  g  be  a  continuous  linear  form,  both  defined  over  36.  The 
problem  to  be  solved  is :  Find  u  €  36  such  that 

(11.1)  V  V  €  36  ,  a(u,v)  =  g(v)  . 

For  the  numerical  resolution  of  this  problem,  we  first  Introduce  a  sequence  of  finite  dimensional 
subspaces  cAl,  c  c  ...  c  <Alj  of  36 ;  then  consider  the  problem  ;  Find  Uj  €  sAtj  such  that 

(11.2) j  V  V  €  <At j  ,  a(Uj,v)=g(v)  . 

The  basic  idea  for  solving  problem  (11.2)  with  a  multigrid  algorithm  consists  In  first 
defining  a  simple  problem  over  the  largest  spaces  and  solving  It,  then  correcting  the  residual 
derived  from  the  solution  of  this  simpler  problem  when  plugged  into  problem  (II. 2)^  by  solving 
problem  (11.2),^  for  lower  values  of  k  <  j.  The  first  step  Is  the  most  important  one  and  relies  on 
the  good  choice  of  continuous,  elliptic,  symmetric,  bilinear  form  b,  called  smoother,  that 
represents  a  in  some  sense  and  is  easier  to  invert.  Let  us  suppose  that  we  have  only  two  grids,  the 
coarse  one  (i.e. ,  )  and  the  fine  one  (i.e.,  </%2  .  ^nd  j  =  2).  The  two-grid  procedure  consists  in 

1  -  m/2  steps  of  smoothing  where  we  solve  m/2  times  a  problem  like  the  following  one  :  Find 

in  cM^2  such  that 

(11.3)  V  V  €  tM2  ,  b(^(p-(p,v)  =  g(v)  -  a(<p,v)  , 

2-  one  step  of  coarse  grid  correction  where  we  solve  only  once  a  problem  like  the  following  one; 
Find?  in  cH,  such  that 

(11.4)  V  V  €  ,  a(?,v)  =  g(v)  -  a((p,v)  , 

and  define  6(p  =  tp  +  ip. 

3-  m/2  steps  of  smoothingas  in  (11.3). 

We  consider  here  just  two  grid  levels,  the  only  reason  is  for  sake  of  simplicity  of  the 
notations,  but  as  in  [B.D.] ,  we  could  consider  the  whole  W-cycles  based  on  more  than  two  levels.  If 
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the  initial  guess  for  the  exact  solution  u  is  u®  ,  after  one  V-cycle  like  the  one  described  previously 
by  the  three  steps  1  -  2-  3-,  the  resulting  solution  Is  u’  and  can  be  expressed  like  a  function  of  u° 
as  follows 

(11.5)  . 

so  that,  after  the  r^^  V-cycles,  the  solution  is 

(11.6)  u'' =  . 

Moreover,  it  is  very  simple  to  note  that  if  u  is  the  exact  solution,  then 
0u  =  u  and  ^u  =  u  , 

so  that,  if  e’’  denotes  the  error  after  the  r^**  V-cycle,  we  derive  from  (11,6)  that 

(11.7) 

Note  that  the  equations  (11.3)  (11.4)  define  affine  operators  0  and  §  but  in  (11.7)  we  can  consider 
these  operators  as  linear  ones  since  they  operate  on  differences  e°  and  e**;  from  now  on,  we  shall 
consider  these  operator  as  linear  ones  while  keeping  the  same  notation. 

As  noted  in  [B.D.],  of  importance  is  the  analysis  of  the  spectrum  of  the  following  eigenvalue 
problem  ;  Find  'i'  in  <H.2  andX  in  [R*  such  that 

(11.8)  V  v  €  cMj  ,  a('l'  ,v)  =  X  b('l'  ,v)  , 

where  b  has  been  scaled  so  that  the  maximum  eigenvalue  X^^^  is  equal  to  1.  Let  us  order  the 
eigenvalues  in  increasing  order  0  <  X^  <  Xj  <  <  Xp  =  X„,^  =  1 ,  (where  P  is  the  dimension  of 

cHj  )  and  choose  relative  eigenfunctions  'P,  ,  'I'g  ..  ..'Pp  .  Of  equal  importance  in  the  analysis  is 
also  the  compatibility  between  the  coarse  space  <Al,  and  the  space  spanned  by  the  first 
eigenvectors . 

More  precisely,  under  the  following  hypothesis 

H  the  space  coincides  with  span  {'1' , ,  '1'2 .  ••  .  (where  p  is  the  dimension  of  <Al,  )  ; 

then  one  can  prove  the  following  theorem 

Theorem  IIJ  :  Under  the  hupothesis  H  the  error  after  the  first  V-cycle  verifies 
a(e’,e’)  <  (1-Xp^,)^"’  a(e®,e®). 

Proof :  First  of  all ,  let  us  recall  that  the  eigenfunctions  ,  n  =  1 ....  ,P ,  form  a  basis  of  cAtj  that 
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Is  orthogonal  for  both  the  forms  a  and  b.  Let  us  span  e°  in  this  basis;  we  get 
=  ; 

then,  due  to  (11.3)  and  (11.8),  we  derive  that  ^''^(e®)  satisfies 
(11.9)  V)  =  ( 1 . 

From  hypothesis  H ,  we  then  derive  that  the  operator0  truncates  the  previous  spectrum  so  that 

e^^V)  =  L„V,(1-X„)"'^^e>„  ; 

then  the  smoothing  procedure  diminishes  once  more  the  spectrum  of  the  error  as  follows 
e’  =  =  Zl,.,  ( 1  -X„)'"  ej  . 

We  deduce  now  from  the  orthogonality  of  the  that 

.«')  =  Z:„"=p„  (I (  e„«)^  e.M-p.l'.)  < 


Remark  II.  1  :  Note  that  the  previous  theorem  is  very  simple  and  is  a  trivial  extension  of  the 
analysis  of  the  multigrid  procedure  in  the  Fourier  space.  Note  also  that  if  the  space  «At,  is  not  so 
well  chosen,  for  instance  if  it  coincides  with  span  {"I'm  -p+ 1>  '^n-p+  2>  multigrid 

procedure  would  not  converge  rapidly  since,  after  the  first  V-cycle,  we  obtain 

and  the  error  remains  important  since  X,  can  be  very  small.  In  fact,  the  method  has  exactly  the 
same  properties  as  the  plain  Jacobi  algorithm. 

It  turns  out  from  the  previous  analysis  that  the  multigrid  procedure,  when  applied  under 
hypothesis  H ,  is  degenerated  since  only  one  V-cycle  is  needed  and  m  is  the  only  important  factor  of 
convergence.  So  in  the  non  trivial  applications  where  the  hypothesis  H  is  not  verified,  we  have  to 
measure  the  actual  situation  between  the  hypothesis  H  and  the  situation  explained  in  remark  II.  1 . 
This  can  be  explained  as  follow  :  the  "rough"  eigenmodes  (  [B.D.)  denotes  this  way  the  with  \ 
close  to  1 )  are  damped  during  the  smoothing  procedure  (  their  components  are  multiplied  by  a 
factor  ( 1  -X„)  )  while  the  "smooth"  modes  (  corresponding  to  small  X„)  remain  almost  constant. 
Under  the  hypothesis  H  these  are  completely  erased  during  the  correction  procedure,  but  if  we  are 
not  in  this  optimal  situation,  they  can  only  be  damped  also  during  this  step. 


-6 


The  general  analysis  proposed  in  [B.D]  allows  for  measuring  the  position  with  respect  to 
both  hypotheses.  Let  us  recall  the  basis  of  their  proof  since  we  shall  extend  it  in  section  V.  To  this 
purpose,  they  introduce  the  various  norms,  defined  for  any  real  index  8,  as  follows 

(11.10)  V4)€cn2.  . 

They  introduce  also  the  function  f 

(11.11)  f(a,b)  =  a® b** (a+b)"^**'*’ =  sup^gjQ (l-x)*x'’  , 

Then  we  derive  that 

(11.12)  ,VT,o<T<e  , 

<  f(m/2,(6-T;)/2)||U|l|,  ; 

let  us  write  now 

III  er''\  III?  =  a(  =  a(  eg""*,  eg""*  -  ?)  =  a(  eg""*,  g""*) 

<  III  eg""*  HI,.,  mg""*  IB,., 

<[sup,....,^iiie*iii,.,/iiie*iii,  1  III  eg""*  III,  III  g""*  III,., , 

SO  that  we  derive 

(11.13)  III  eg""*  III,  <  [sup,,  ^  III  e*  III,.,/  in  e*  ni,  i  in  g""*  111,.* . 

which  is  valid  for  any  &  >  0;  we  deduce  from  (11.12)  that 
111  III,  <  f(m/2.6/2)  til iii,.^ 

<  /(m/2.^/2)  III  III,  [sup^g  ^  Jll  (5i|>  III, .p/ III  3«|>  III,  ] 

<  /(m/2,&/2)  (sup^,  ^  ^  III  e«|,  |||,.p/  III  enp  III,  f  III  |||„p 

<  Km  ,&)  [sup,^^  ^  ^  III  Cip  |||,_p/  III  e«p  III,  f  III  (p  III,  . 

Defining  </!,•*■  as  follows 

oM,'*'  =  { Ip  €  <^2  ,  such  that  V  ip  €«At,  ,  a(ip,ip)  =  0  }  , 

( that  coincides  with  the  range  of  Q  )  and  by  minimizing  the  right-hand  side  over  they  state  the 
following 

Theorem  II.  1  :  Assume  there  exists  a  constant  k  ^  1  and  (x  >  0  such  that  for  any  ip  € 

III  *  III,..  <K*"  111*111,  , 
then 

(II.  M)  a(e\e’)  ^  •y  a(e®,e°)  , 
where 
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(11.15)  ^  =  [(k-O/k]^"’  if  m<o((K-1) 

■y  =  (K*r(m,a))^  if  m>oc(K-l). 

11.2  Formulation  of  the  spectral  element  discretization. 

Let  us  turn  now  to  the  position  of  the  problem.  We  consider  here  the  simple  test  problem 
over  the  interval  A  =  ]  - 1 , 1  [ ;  Find  u  such  that 

(11.16)  -  ^  ■  over  A 

provided  with  homogeneous  Dirichlet  boundary  conditions 

(11.17)  u(-l)  =  u(1)  =  0. 

where  f  is  a  given  force.  This  problem  is  very  simple,  but  it  allows  a  statement  of  the  basic 
features  of  the  multigrid  algorithm  and  an  understanding  of  why  the  method  works. 

The  spectra)  element  method  for  approximating  the  solution  of  (11.16)  consists  in 
discretizing  the  space  of  acceptable  functions  by  a  subspace  of  piecewise  polynomials.  More 
precisely,  given  a  pair  h  =  (K,N),  we  first  break  the  interval  A  into  K  disjoint  subintervals  of 
comparable  sizes 

^  -  *^ka1  >  ^k  =  1  ^  ^  *^k  t  ' 

then,  we  choose  for  space  of  approximate  functions  a  subspace  xjj  of  Hq(A)  consisting  of  all 
piecewise  polynomials  of  degree  <  N, 

(II.  18a)  xJ  =  Yl'nH^(A)  , 
where 

( 1 1 . 1 8b)  yJJ'  =  { (p  such  that  €  IPn(  A^^) ) 

and  IPn(A,^)  denotes  the  space  of  all  polynomials  of  degree  <  N  on  A^ .  We  remind  that  contrarily  to 
the  finite  element  approximation,  the  convergence  is  achieved  by  increasing  the  degree  of  the 
polynomials  N  and  not  refining  the  mesh. 

The  discrete  problem  starts  from  the  variational  formulation  of  problem  (II.  16)(ll.  1 7) 
that  is :  Find  u  in  Hq( A)  such  that 
(11.19)  Vv€H^(A).  a(u,v)  =  (f,v)  . 

where  (.,.)  denotes  the  L^(A)-scalar  product  and  a(.,.)  denotes  the  following  bilinear  form  defined 
over  Hq(A)  as  follows 
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V  (p,tv  €  Hg(A)  ,  a((p,t|/)  =  . 

Then  we  construct  the  numerical  scheme  by  discretizing  with  Gauss-Lobatto  quadrature  formulas 
the  various  integrals  present  in  (II.  1 9)  and  restricting  the  space  of  test  functions  to  xjl . 

This  results  in  the  following  problem  ;  Find  u^,  €  ,  such  that 

(11.20)  Vv^tx;  ,  ,v„)  =  (f,v,)“ot  . 

Where  the  discrete  forms  are  defined  as  follows 

V  (P.4.  €  y;;  ,  =  21^,1  ^/2  lL  oH  . 

.  ah.oL^'P-'l')  =  (41,  .»l',)h.GL  • 

Here,  the  ,  and  the  are  the  weights  and  nodes  of  the  Gauss  Lobatto-Legendre  formula  with 
N+ 1  points  and  the  collocation  points  t;nk  defined  by  =  a,^  +  +  1  )b^^/2  .  W  recall  here 

that  the  integration  formula  is  exact  on  IP2n_i  that,  contrary  to  the  appearances,  a||'  does  not 
depend  on  h  nor  N  since  =  a  (  see  e.g.  [M.P.]). 

The  algebraic  system  that  has  to  be  solved  is  derived  by  choosing  the  values  of  the  unknown 
function  u^  on  the  collocation  points  and  representing  u^  in  the  basis  of  the  interpolant  basis  h,^„ 
defined  by  ;  h,^  „  is  the  only  element  of  yJJ'  such  that 

(11.21)  V  \/!;,„)  =  6„6„„  . 

The  matrix  system  that  has  to  be  solved  can  be  written  as  follows:  Find  u  =  (u||''’)  such  that 

(11.22)  Au  =  g  , 

where  A  is  the  stiffness  matrix  with  entries  equal  to 

(11.23)  (2/b^)L;".,  Zlo  dhj„/(Jx)  (C.,^)  , 

with  =  P„b^/2  and  21  denotes  the  direct  "stiffness  summation"  ,  while  g  is  related  to  the 
forcing  term.  We  refer  to  [R.P.]  so  as  to  [M.P.]  for  more  details  on  the  derivation  of  this  matrix. 

This  numerical  method  is  proved  to  work  ve^'y  well  in  a  great  number  of  interesting 
problems  that  include,  for  instance,  the  full  Navier-Stokes  problem  (see  for  instance  [M.P.], 
[M.P.R.  1  ],  [M.P.R.2]  or  [R])  and  is  numericaly  competitive  and  implementable  on  a  parallel 
medium  grain  paradigm  (see  for  instance  [F.R.D.P.]). 


Let  us  turn  now  to  the  multigrid  algorithm  for  solving  iteratively  problem  (11.22).  As 


explained  in  [R.P.]  ,  nested  spaces  are  related  to  spaces  of  polynomials  with  lower  degree  (say 
M  =  N/2)  and  the  smoother  is  simply  the  Jacobi  preconditioner  B  that  is  proportional  to  the 
diagonal  part  of  the  matrix  A  and  normalized  in  such  a  way  that  the  highest  eigenvalue  of  B'^A  is  1 . 
We  shall  analyze  in  the  two  following  sections  the  properties  of  this  preconditioner  and  explain 
why  the  numerical  method  works  so  well  as  presented  in  [R.P.] . 

In  this  paper,  we  shall  use  two  grids  only  though  the  analysis  can  be  performed  with  no 
extra  difficulty  than  the  one  of  comprehension,  and  we  shall  use 
cM,  =  and  ^^2  =  ^11  • 

The  problem  enters  in  the  general  theory  of  [B.D.]  since  the  numerical  problem  involves  bilinear 
forms  and  matrices  that  do  not  depend  on  h. 
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III.  Analusis  of  the  convergence, of  the  multiqrld  algorithm  In  the 
case  of  one  element. 

As  announced  in  the  title,  we  shall  assume  for  the  moment  that  the  discretization  is  applied 
on  only  one  domain  to  problem  (11.16).  For  sake  of  simplicity,  we  shall  drop  out  the  second  index 
that  corresponds  to  the  subelement  characterization  in  the  notations.  The  situation  is  simple  here 
as  soon  as  we  have  explained  the  properties  of  the  diagonal  matrix,  but  this  simple  example 
permits  enhancement  of  the  main  features  that  allow  for  a  rapid  convergence  of  the  algorithm.  The 
problem  can  be  written  as  a  pure  collocation  scheme  as  already  noted  in  [M.P.].  Indeed,  by  taking 
'^h  =  for  n  =  1 ...  ,N- 1 ,  in  (11.20)  and  using  twice  the  exactness  of  the  Gauss-Lobatto  formula, 
we  derive 

=  '  Ia  =  -  (  U^*  ,  h,),^^  =  -  i 

besides  we  note  that 

so  that  the  problem  actually  verified  by  u^  is :  Find  u,,  in  xj^  such  that 
V  n,n=l,..,N-l  ,  -u;'(g  =  f(g  . 

This  consists  in  a  pure  collocation  procedure  to  solve  the  initial  problem. 

In  order  to  analyze  the  multigrid  algorithm ,  let  us  first  compute  exactly  the  diagonal  part  of 
the  stiffness  matrix;  that  is  here  simply  to  the  matrix  with  entries  equal  to 
=  .fori.jrl . N-l  , 

^Oj  ■  ^Nj  ~ 

Due  to  the  exactness  of  the  Gauss  Lobatto  formula,  we  easily  obtain  that,  for  i  =  1 . N- 1 , 

Aji  =  [hj'(x)]^  dx  =  -  h."(x)  hi(x)  dx  , 

and  here  again,  by  using  the  exactness  of  the  Gauss  Lobatto  formula,  we  derive  from  (11.2 1 ) 

As  already  derived  in  [GHO,  formula  (7.^)] 

Lemma  111.1  :  The  diagonal  of  the  stiffness  matrix  verifies 
(III.1)  Aii  =  p|N(N+l)/[3(l-(^[')^)J. 

Proof  :  Let  us  drop  out  in  this  proof,  the  superscript  N.  First  of  all,  we  note  that,  from  the 
definition  (11.21)  of  h,. ,  we  have 
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hjCx)  =  _ _  =  ocj  nfrj  (x-Cj)  (1-x^) 

n,;,  (ti-yccf-i)  "■ 

where  Wj  is  a  non  zero  constant.  It  is  well  known  (see,  e.g  [D.R.])  that  the  internal  nodes  of  the 
Gauss  Lobatto  fomula  verify 

(111. 2)  Ln(x)  =  c,  n-:I(x-t;p  . 

where  is  a  non  zero  constant.  Hence,  we  can  write 

(x-J;;)  hjfx)  r  cx,  ( 1  -x^)  L^(x)  =  oc.  ( 1-x^)  Lj;,(x), 
and  after  taking  the  derivative  of  both  sides,  we  obtain 
(x-f:;)  h;(x)  +  hj(x)  =  otj  (d/dx)[(  1  -x^)  ^(x)]  . 

Let  us  recall  now  the  eigenfunctions  property  verified  by  the  Legendre  polynomials  (see,  e.g. 
[D.R.]) 

(111. 3)  (d/dx)[(l-x^)  L,;,(x)]  = -(N)(N  +  1)  Ln(x)  ; 
we  derive  that 

(111.4)  (x-{;.)  h:(x)  +  hi(x)  =  -  oCjN  (N  +  1)  Ln(x)  ; 
plugging  now  x  =  i;.  in  this  equation  yields 

(111. 5)  1  =-0CiN(N+l)LN(g  . 

Besides,  by  taking  the  derivative  of  (III. 4),  we  obtain 

(111. 6)  (x-?;i)h;'(x) +  2h:(x)  =  -a.N(N+1)L;,(x)  , 
plugging  also  here  x  =  ^.  in  this  equality  and  using  (111.2),  we  get 

(111. 7)  h:(t;i)  =  0. 

Let  us  multiply  now  (111,6)  by  ( 1  -x^)  and  take  the  derivative  of  the  resulting  equation,  we  deduce 
(x-g  (d/dx)[(l-x2)  h;'(x)]  +  3(  1 -x^)  h."(x)  -  4xh;(x)  =  oc,  (N  +  1)^  Ln(x)  . 

Finally,  plugging  one  more  time  x  =  ,  we  derive  from  (III. 5) 

3(i-^;f)h;'(g  =  -N(N+i)  , 
and  the  lemma  follows  from  the  value  of  Ajj  . 

Remark  III.1  :  Note  that,  as  a  consequence  of  (1 1 1.7),  we  have  proved  here  that  the  Lagrangian  at 


point  i;.  has  its  maximum  at  {;, . 
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Assoclated  with  the  matrix  A  is  the  bilinear  form  a,  and  in  the  same  way,  associated  with  B , 
the  normalized  diagonal  of  A  can  be  defined  a  bilinear  form  bfj  smoother  of  the 

multigrid  algorithm.  Lemma  III.  1  proves  that  the  smoother  we  have  introduced  is  proportional  to 
the  bilinear  form  b|J|  ql  defined  for  any  ip  and  41  in  xjj  (which  here  is  simply  IP,^(A)  H  Hq(A)) 

(1 1 1.8)  b,%L(<p.<i/)  =Ll.o  Pi<p(^i)  4'(«;i)  d-tfr’  • 

It  is  interesting  to  note  that,  due  to  the  exactness  of  the  Gauss-Lobatto  formula,  we  can  rewrite 
GL  ®  continuous  way  since  for  any  (p  and  <p  in  xJJ  ,  the  function  [«p  (p  ( 1  ]  is  still  a 

polynomial  and  belongs  to  1P2n.2(A),  so  that  we  have 

(111. 9)  =  6('P.4')  =  Ia  'p(x)  «p(x)  ( 1 -x^)"’ dx  . 

Let  us  now  analyze  the  eigenvalue  problem  (11.8)  or  more  precisely  the  eigenvalue  problem 
associated  to  b.  The  situation  is  here  very  simple  since  we  have  an  exact  expression  for  the 
solution  to  this  problem 

Lemma  111.2  ;  Let  us  define  for  eng  integer  n  ,  1  <  n  <  N- 1 

(111.10)  'l'„(x)  =  ( 1 -x^)  L;(x)  ; 
then  we  have 

(111.11)  Vv^xJ  ,  a('P„,v)  =  n(n*l)B('4/„,v)  . 

Proof  :  Let  us  first  remind  the  following  standard  formula  on  the  Legendre  polynomial  (see,  e.g. 
[D.R.,Chap2,§7]) 

(111.12)  Vn€lN,  ((1-x^)l;(x))' +n(n+l)L„(x)  =  0. 

Let  us  compute,  for  any  v  in  x|J  , 

a(^^„  ,v)  =  'l';(x)  v'(x)  dx  =  (( 1  -x^)  l;(x))'  v'(x)  dx  ; 

using  (III.  1 2)  and  integrating  by  parts,  we  obtain 

a('l'n  ,v)  =  -  n(n+ 1 )  L„(x)  v'(x)  dx  =  n(n+ 1 )  L' (x)  v(x)  dx 

=  n(n+ 1 )  ( 1  -x^)  L;(x)  v(x)  ( 1  -x^)"’  dx  =  n(n+ 1 )  b('l'„ ,  v) ; 

this  ends  the  proof  of  Lemma  III.2  . 

It  is  important  to  note  that  in  this  simple  example  the  eigenvalues  are  well  known  and, 
moreover,  that  the  first  M  of  them  span  exactly  the  space  Pm^/A)  fl  Hq(A).  As  a  first 
consequence,  we  can  state  that  the  normalized  form  that  will  be  used  as  a  smoother  is  defined  as 
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follows 

(111.13)  Vtp,ii/€Xj  ,  b((p.4/)  =  N(N-l)b(0,.j/)  =  N(N-1)J^  <p(x)  v|/(x)/(l-x^)  dx  . 

We  demonstrate  the  case  where  the  simple  hypothesis  H  of  section  1 1  is  valid,  and  we  can  state  now 
Theorem  111.1  :  Let  u°  denote  the  initial  guess  in  the  multigrid  algorithm  applied  to  problem 
(11.20)  in  the  case  of  one  element,  and  u’  denotes  the  solution  obtained  in  the  after  m  smoothing 
and  one  correction.  This  solution  converges  to  the  exact  solution  u  as  follows. 

a(u-u\  u-u’)  =  [1  -  (N+2)/-4(N-l)]^'"a(u-u°,  u-u®)  . 

Proof  ;  This  is  a  simple  corollary  of  Theorem  11.1  and  lemma  111.2  since  the  first  (N/2)-1 
eigenvectors  span  the  coarse  space  ,  it  follows  that  the  eigenvalues  of  problem 

(11.8)  are  n(n+ 1  )/(N- 1  )N. 

Remark  111.2  :  Note  here  that  the  correction  on  the  coarse  grid  needs  only  be  done  once  and  that 
there  is  no  optimal  choice  for  the  number  of  smoothings  since  the  convergence  is  proportional  to 
this  number.  This  is  actually  in  accordance  with  the  numerical  simulation  of  [R.P.]  as  appears  in 
table  1  of  that  paper. 

Remark  111,3  :  Let  us  point  out  the  fundamental  reasons  that  give  this  rapidity  to  the  algorithm. 
They  are  hidden  here  due  to  the  simplicity  of  the  eigenvalue  problem.  First  of  all,  even  if  this  is 
not  of  major  importance,  the  matrix  B  is  a  good  preconditioner  of  the  matrix  A.  Indeed,  the 
condition  number  of  B'’a  is  order  as  opposed  to  the  condition  number  of  A  which  is  order  N^. 
This  will  be  also  the  case  for  other  problems  and  has  already  been  noticed  by  (H.]  in  a  different 
context  for  an  application  to  conjugate  gradient  algorithms.  This  can  be  viewed  as  an  inverse 
inequality  or  a  Jackson  type  one  since  the  weighted  L^-type  norm  associated  to  the  bilinear  form  B 
is  com  pared  to  the  H  \  A )  -  nor  m  as  fol  lows 
(lll.M)  V  (p  €  IPn(A)  DH^CA)  .  ||«p||?^^  <N(N-l)b(<(),ip)  . 

Secondly,  the  other  property  that  is  very  important  in  the  multigrid  algorithm  is  that  the  factor  k 
as  defined  in  Theorem  11.1  is  also  bounded  by  two  constants  independent  of  N.  Indeed,  we  can  state 
V  ip  €  xJJ  such  that  a{J|  g^^tp  ,v)  =  0  for  any  v  in  xJJ^^ 
we  have 

(III.  15)  b(<p,(p)  <  a((p,(p)  , 
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where  the  constant  k  is  bounded  by  4(N- 1  )/(N+2)  as  can  be  derived  from  lemma  1 1 1.2  . 

Let  us  generalize  now  the  results  obtained  we  have  obtained  in  this  very  simple  situation  to 
the  case  of  a  multielement  discretization. 


We  begin  here  also  by  analyzing  the  properties  of  the  diagonal  of  the  stiffness  matrix  A.  We 
immediately  note  that  there  are  two  kinds  of  diagonal  elements  in  this  matrix  ;  those  that 
correspond  to  internal  points,  i.e.  that  involve  the  scalar  product 
(IV.  1  )  Aj,  =  1  ^nsO  9k,n  [dhj  j/dx  dhj|/dx](tnj^) 

with  i  =  1  ,  and  those  that  correspond  to  interface  points,  i.e.,  that  involve  the  scalar 

product  (IV.  1 )  for  i  =  0  or  N.  The  first  kind  of  diagonal  elements  are  the  same  as  those  involved  in 
the  previous  section.  Indeed,  the  corresponding  Lagrangian  interpolants  vanish  at  the  interfaces 
and  also  on  any  subinterval  that  does  not  contain  the  point  ;  therefore,  from  Lemma  III.  1 ,  we 
can  state  that 

V  i,i=  1,..,N-1  .  VC  ,l>  =  1,..,K  .  A„=(2/b,)p.N(N+l)/[3(1-(j;i)^)l  . 
or  again,  thanks  to  a  simple  change  of  variable 

(IV.2)  Vi,i=  1...,N-1  .  VC  ,C  =  1,...K,  Aj,  =  p,,  N(N+ D/CSfJ;,*- a,)(a,^,-Li,,)]  . 

For  the  interface  terms  we  have 
Lemma  IV.  1  :  For  i  =  0  and  any  C  =  2,..,K  we  have 
(IV.3)  A.,=  t(b,.ir’+(b,r’](N^N  +  l)/3  . 

Proof ;  As  already  used,  the  exactness  of  the  Gauss-Lobatto  formula  gives,  for  i  =  0  and  C  =  2 ,..  ,K 
(  or  i  =  N  andC  =  1 ,..  ,N- 1 ,  due  to  the  direct  stiffness  summation) 

(IV.4)  A.,  =  dl  *  [dho/dxf  (0  dL  . 

A  simple  change  of  variables  and  the  use  of  the  symmetry  of  h^  and  h^  yields 
(IV.5)  =  [(2/b,.,)+(2/b,)l  [dh„/dxl*(0  dt,  . 

Let  us  compute  the  integral  on  the  right-hand  side  of  this  equation.  First  we  have 


hN(x)  = 


(x+1)  (x-?;,) 


From  (III.  1 3)  we  then  get  that 

hN(x) = (x+1)l;(x)/2L;(1)  , 
so  that,  after  integration  by  parts,  we  obtain 
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1a  dx  =  -  (2  L,;,(l))-^  [J^  (x+1)  l;(x)  (d^/dx^)((x  +  1)  L;;(x))  dx 

-2L;,(1)  (d/dx)((x+l)  Ln(x))(1)  ]  ; 

here  again,  the  use  of  the  exactness  of  the  Gauss-Lobatto  formula  to  compute  the  integral  on  the 
right-hand  side  yields 

( IV,6)  [dh/dx]2(x)  dx  =  -  (2  L,;,(  1  ))•’  [(d^/dx^)((x+ 1 )  l;(x))(  1 ) 

-  (d/dx)((x+l)  L,;,(x))(l)]  . 

It  is  an  easy  matter  to  note  that 

(d/dx)((x+ 1 )  Ln(x))  =  (x+  1 )  Ln  (x)  +  L;(x)  . 

(d^/dx^)((x+ 1 )  L^Cx))  =  (x+ 1 )  Lj;,"  (x)  +  2  Lj;j  (x)  , 
from  (111.12)  writen  in  the  form 

( 1  -x^)  L,;;  (x)  -  2xL;,(x)  +  N(N+  l )  L^Cx)  =  O  ; 
we  derive  that 

L;,(1)  =  N(N+l)/2  . 

L,;;(1)  =  (N-l)N(N+l)(N+2)/8  . 

L,;,"(l)  =  (N-2)(N-l)N(N+l)(N+2)(N+3)/48  ; 
this  gives 

(d/dx)((x+l)  Ln(x))(1)  =  N^(N+1)V4  , 

(d^/dx^)((x+l)  l;(x))(1)  =  (N-l)N^(N+l)^(N+2)/24  . 

Plugging  this  in  (IV. 6)  and  using  the  relation  (see  [D.R.]) 

(IV.7)  Pn  =  2/(N+1)N  , 
we  derive 

(IV. 8)  [dh,^/dx]^(x)  dx  =  (N^+  N  +  1  )/6. 

The  lemma  follows  then  from  (IV.5). 

From  (IV. 2)  (IV. 3)  and  (IV, 7),  we  derive  that  the  bilinear  form  b|||gL  that  defines  the 
smoother  and  is  associated  to  the  normalized  diagonal  part  of  A  is  proportional  to  the  bilinear  form 
^h6L  defined  for  any  ip  andip  inxjj'  as  follows  (remind  (1 1 1.8)) 

(IV, 9)  bh  QL(ip,<l»)  =  List  B,(ip,(i()  , 
where 
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+  Po.i'P^V*)  N  +  0/  (b,)*]  +  Pn.i'P^W  '•'(W  t(N^  N  +  0/  (b,)2] } . 

rN 

In  the  two  next  lemmas,  we  shall  analyze  the  eigenvalue  problem  between  a  and  b|,  that 

N  N 

will  allow  first  to  estimate  the  normalization  factor  between  b^  ol- 

Lemma  IV.2  :  For  any  <f  in  Xj  ,  we  have 

(IV.  10)  J^<p'^(x)dx<(4N^+2N2+3N-1)N/3(N^+N+1)Bh5L('P.'l')  • 

Proof ;  To  any  tp  in  xJJ  let  us  associate  the  element  (p^  defined  as  follows 
(IV.  11  )  <Po  =  <P  -  Lf,2  [%.«  +  • 

so  that  (Pq  and  (p  coincide  on  the  interior  collocation  points.  Using  the  inequality 
(a+b)^  <  ( 1  +a"’)  a^  +  ( 1  +a)  b^  , 
we  deduce  from  (IV.  11), 

(IV.  12)  dx  <(!+«"’)  1a, dx  +  (1+a)  [<p(a,)ho,  +  <P(a,^,)hN,,]^(x)  dx. 

It  is  an  easy  matter  to  note  that  the  restriction  of  ip^  to  any  A,  belongs  to  ^^(A,)  H  Hq(A,)  so 
that  the  lemma  1 1 1.2  and  a  simple  change  of  variables  yields 

1^,  dx  <  N(N-1)  Zu!  Q;,t  a,)(a*+brMr^ 

<  N(N-l)  ^1*^’  Pi,!  a»)(at+b,-Ci,,)]’’. 

Besides,  from  (IV.8),  we  derive  that 

1a,  thi,,!*  <*)  *>  =  1a,  ["«.,]*('')  *<  =  (N**  N  .  D/Jb,  . 
as  following  the  same  lines  as  in  the  proof  of  lemma  IV.  1 ;  we  get,  for  any  N  >  2 
1ai  bo,<  (x)  dx  <  1  /3b,  . 

F rom  ( IV.  1 2)  we  then  deduce 

f‘(x)  dx  <  ( I  *«-')N<N- 1 )  L.".','  ?i,  ♦“(ti,)!  (K;r 

♦  (!*«)((N®tN.2)/Jb,ltl»(a,)|*  t|*(a,.,)l*). 
choosing  now  a  =  3(N-1)(N*+N*  1)/(N+1)(N**N+2);  it  follows  from  (IV.7)  and  (IV.9)  that 
1^  «’^(x)  dx  <  (4N^2N^3N-1)N/3(N*+N+1)  . 

and  the  lemma  follows. 

Remark  IV.  1  :  The  estimate  in  the  previous  lemma  provides  a  less  precise  characterization  of  the 
eigenvalue  problem  than  the  one  we  could  get  in  lemma  ill.1 ,  but  we  note  that  the  highest 


eigenvalue  involves  the  same  asymptotic  order  as  in  the  previous  section  and  this  will  be  enough 
for  our  purpose.  The  important  fact  is  that  the  result  is  Independent  of  K  and  of  the  ratios  between 
the  various  subinterval  length  b, .  As  we  shall  see  in  what  follows,  this  will  result  in  a  multigrid 
algorithm  that  will  work  as  well  in  any  case  of  number  of  subelements.  Besides,  note  that  the 
smallest  eigenvalue  of  problem  (11.8)  scales  like  (independently  of  N)  such  that  the  condition 
number  of  B"’a  behaves  like  (KN)^  ,  in  accordance  with  the  finite  element  equivalent  (when  N  is 
order  1 )  and  proves  that  the  conjugate  gradient  algorithm ,  when  preconditioned  by  B .  has  a  rate  of 
convergence  which  behaves  like  1  -  c/KN  and  depends  on  both  K  and  N!!!  This  is  of  importance  when 
we  compare  the  preconditioned  conjugate  gradient  with  the  multigrid  algorithm.  We  refer  to  the 
thesis  of  E.  M.  RONQUIST  [R]  for  numerical  evidences. 

It  is  an  easy  matter  to  derive  from  lemma  IV. 2  that  the  normalized  form  b^  defined 
from  bf,  QL  by  multiplication  by  a  factor  of  order  (4/5)N  .  The  other  property  that  is  important 
for  the  analysis  of  the  multigrid  algorithm  deals  with  the  space  of  those  elements  (p  of  X{J 
that  verify 

(IV.  13)  V^€X;"2  .ah%L(«P.<k)  =  0  . 

Lemma  LV.3  :  For  any  «p  in  we  have 

(IV.M)  bJj5L('P.'P)  <4(4N^2N^3N-l)/3(N^N+l)(N42)a(«,q))  . 

Proof :  Let  (p  belong  to  ,  and  let  us  define  for  ?  =  2,..  ,K  the  element  ip,  of  x{J'*  by 
g/,€X^’  .  Vk=  I,...,K>1,  <pA)  =  8|k  • 

It  is  an  easy  matter  to  compute  that 

,  for  X  in  A,_,  , 

<pj(x)=  -I/b, ,  for  X  in  A,  , 

0  ,  forx  inA^  ,k?s^  andk;^?-!  . 

Using  this  function  in  ( IV.  1 3) ,  we  derive  that 

V«=2,...K,  (1/b,.,)(A)- A-i»^(’/M(«(a,)-«(a,^,))  =  0  j 
recalling  now  that  (p(ao)  =  =  0,  we  deduce  that  in  fact  41  vanishes  at  any  interface  so  that 

(p,^^  is  an  element  of  D  Hq(A,).  The  use  of  (111.15)  over  each  A,  proves  that 
B"6l(<p,«p)  =  LJ.,  BiCf.*)  ^  4/N(N+2)  a(«.«)  . 
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and  the  lemma  follows  by  recalling  the  normalization  factor  of  b^  ql- 

As  a  simple  consequence  of  the  previous  lemmas  and  (11.1  1)(II.14)  we  derive  that  the 
multigrid  algorithm  converges  toward  the  numerical  solution  with  a  speed  independent  of  both  N 
and  K ,  indeed,  we  have 

Theorem  IV.  1  :  The  multigrid  algorithm  based  on  the  Jacobi  preconditioner  converges,  and  at  any 
V-cycle  with  m  smoothing  as  detailed  in  section  II  ,  the  rate  of  convergence  is  given  by 
|;^ef(x)dx  <  Y  j^eQ^(x)dx  . 

where 

U  =  [(k-D/kI^'"  ,  for  m  <  (k-1)  . 
and 

V  =  (Kf(m.l))^  ,  for  m>(K-l), 

with 

K  =  4(4N^  +  2N^  +  3N-l)/3(N^+N  +  l)(N+2)  =  ( 1 6/3)  ( 1  +0(N’’))  . 

Remark  IV. 2  :  Let  us  note  that  the  convergence  rate  that  we  have  theoretically  obtained  is 
independent  of  K  and  the  sizes  of  the  subintervals  and  does  not  deteriorate  when  N  increases;  this  is 
in  accordance  with  what  is  numerically  ob^rved  in  part  one  of  this  paper  [R.P.J.  However,  the 
rather  rough  estimate  we  gave  for  the  highest  eigenvalue  in  (lY.  10)  provides  a  rather  too  high 
estimate  for  the  convergence  parameter  (  close  to  0.8 1  when  N  is  large  enough  ).  Note  also  that  the 
optimal  choice  of  parameter  4  <  m  <  5  is  close  to  what  can  be  observed  numerically.  As  in  the 
analysis  of  the  multigrid  algorithm  when  a  finite  element  method  is  used,  this  optimal  parameter 
is  lower  than  the  actual  one.  By  using  negative  order  of  the  norms  defined  from  b  and  a,  we  could  fit 
more  closely  to  the  experiments  for  this  last  result.  This  would  be  obtained,  however,  at  the  price 
of  a  much  more  technical  proof  and  would  not  really  be  worthy  since  the  main  conclusion  is  the 
independence  of  the  convergence  rate  with  respect  to  the  parewneters  of  the  discretization. 


The  previous  chapter  was  devoted  to  the  analysis  of  the  multigrid  algorithm  when  applied  to 
the  very  simple  equation  -u^^  =  f.  Of  interest  of  course  is  to  know  that  the  same  conclusions  hold 
true  in  more  complicated  situations.  We  shall  extend  here  the  analysis  to  the  case  of  the  equation 
(V.l)  -(au,)^  =  f  . 

where  a  is  a  function  of  x  such  that  there  exists  two  constants  a"  and  oc*  with 
(V.2)  Vx€A,  0<oc'<  a(x)  <  a*  , 

and  also  such  that  a  is  in  the  Sobolev  space  A).  We  shall  assume  here  that  the  domain  is  not 
decomposed  into  subdomains,  leaving  this  analysis  to  a  forthcoming  paper  as  well  as  the  analysis  of 
the  multidimensional  case.  The  analysis  provided  here  is  inspired  from  the  reference  [B.V.]. 

It  is  standard  to  note  here  that  problem  (V.  I )  can  also  be  stated  in  a  variational  formulation 
like  (II.  19)  with  now  a  defined  as  follows 
(V.3)  V  €  HqCA)  ,  a((p.t|i)  =  f^0((x)  (p*(x)  ♦‘(x)  dx  . 

It  is  rather  well  known  also  that  the  general  spectral(element)  discretization  of  the  equation 
consists  in  first  choosing  a  discretization  parameter  n  €  W*,  then  :  Find  u„  €  X"  =  IP^CA)  D  H^CA) 
such  that 

(V.4)  Vv„€X"  ,a'’(u„.v„)  =  (f.v„)„  . 
where  the  discrete  form  a"  is  defined  here  by 
(V.5)  Vip.^reX",  a"(<p.»i»)  =  (a  Ip,  . 

The  nested  spaces  for  the  multigrid  algorithm  are  exactly  the  same  as  in  the  previous  sections,  i.e. , 
and  cMj  =  X^  ,  and  the  strategy  also  based  on  the  use  of  the  Jacobi  preconditioner  as  a 
smoother.  The  corresponding  bilinear  form  b**  is  deduced  from  the  following  6*^  after  normalization 
of  the  maximum  eigenvalue  of  problem  (11.7):  S'*  is  defined  by 

(V.6)  v«.^€xv  6^ <(»,♦)  =  pH  - 

where  corresponds  to  the  diagonal  entry  of  the  stiffness  matrix  A,  equal  to 

and  we  recall  that  h”  is  the  Lagrangian  interpolant  at  the  point  .From  hypothesis  (V.2)  and  the 
exactness  of  the  Oauss-Lobatto  quadrature  formula,  it  is  simple  to  derive  that 
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(V.7)  V  (p  €  ,  a'^(ip,i(i)<  a*  aCip.ip)  <  (X*  N(N-1)  B((p,(p) 

<(3aVa-)(N-l)/(N+1)  6'^((p,(p). 

where  we  recall  that  bCip.ip)  =  <p(x)  4)(x)/(l-x^)  dx.  As  a  result,  we  derive  that  the 

normalization  of  satifies 

(V.8)  b^  =  (3aVa')(N-l)/(N  +  1)  6^. 

From  the  eigenvalue  problem  :  Find  in  andX  in  IR*  such  that 
(V.9)  V  V  €  X\  a^^.v)  =  X  b^^'.v) 

that  possess  N- 1  eigenvalues  0  <  X^  ^  X^  <  ...<  X^_,  <  1 ,  the  eigenvectors  of  which  are  chosen 
normalized  with  respect  to  the  norm  derived  from  b^  we  define,  as  in  (II.  10)  the  ||| ,  IHg^  norms  for 
any  0  €  IP  by 

(V.  1  0)  V  V  €  XV  III  V  llll^  =  LH'Ji  x; 

Here  the  bilinear  forms  associated  to  the  definition  of  the  system  and  to  the  smoothing  depend  on  N. 

The  definitions  of  the  smoothing  operator^  and  the  correction  operator  0^  have  to  be  precised.  Let 

us  do  this  in  the  simple  case  where  only  two  grids  are  used.  Instead  of  (11.3),  the  smoothing 

procedure  consists  in  :  Find  in  X^  such  that 

(V.ll)  Vv€XV  bV^%-4),v)  =g^v)-a^(p,v)  . 

while  the  correction  procedure  consists  in :  Find  tp  in  X  such  that 

(V.  1 2)  V  V  €  =  g^v)-a'*(,p.v)  , 

and  define  0^(p  =  (p  +  Then  as  explained  in  the  general  case,  during  a  V-cycle  with  m/2 
smoothings  down  and  up  the  error  e°  is  changed  in  e’  as  follows 
(V.13)  e’  =  (^)'"'2($N(0N)m/2gO 

Before  entering  in  the  details  of  the  analysis  of  the  decay  rate  of  the  error,  let  us  state  some 
results  of  general  interest.  First  of,  let  x  be  in  IR°(A)  be  given,  r(x)  and  r,^(x)  in  (P^fA)  be  such 
that 

(V.M)  V  «j/ €  IPj(A)  ,  b(r(x),'p)  =  b'*(r„(x),<p)  =  a''(x,'l»)  . 

where  b  is  defined  in  (III.  1 3).  It  is  a  simple  consequence  of  (V.9)  to  derive  that 

III  ''n^x)  III?.n  =  b"(rN(x),rN(x))  =  a^x-r^Cx))  <  III  X  III2,n  111  IIIo,n  - 


whence 


Besides,  from  (V.7)  and  (V.  1 4),  we  have 

b(r(y).r(x))  =  b^^N(x)  .r(y))  <  ^^(x))  b^^(x)  ,r(x))]’^^ 

<  (a*Vo(-)’^2  [bN(r^(x),rN(x))  b(r(x)  .r(x))]’'^  ; 
finally,  we  derive  that  the  solution  r(x)  of  (V,14)  has  the  following  stability  property 

(V.15)  b(r(x),r(x))  <c|||  xlll2,N  • 

Then  let  us  state  some  approximation  result  the  proof  of  which  will  be  presented  in  the  appendix. 
Lemma  V.l  ;  Let  p  defined  over  A  and  such  that 
(V.  16)  0^(x)/(l-x^)  dx  <  oo  , 

and  n  be  a  positive  integer.  The  solutions  (P(q)  €  Hq(A),  tp„(g)  €  IP®(A)  defined  by 
(V.l  7)  Vv€H’(A),  a(4)(Q),v)  =  b(Q.v)  , 

(V.l  8)  Vv€(PJ(A)  ,  a'’((p„(p).v)  =  B(p,v)  ; 
then  the  following  approximation  results  holds 
(V.l 9)  I  <p(q)  -  ((I„(q)  I?  <cn"^b(p,p)  . 

(V.20)  B((p(p)  -  (p„(p),(p(p)  -  (p„(p))  <cn*^6(p,p)  . 

Let  us  denote  by  e  the  term  (^)"’^^e®  ;  derive  now  as  in  (II.  1 2)  that, 

(V.21)  ll|e’|||,,N<f(m/2,1/2)l||(?'e||lo.„  , 

(V.22)  |||c|ll2.N<^nn/2,l/2)|||e®|||,N  . 

From  the  definition  of  (3*^,  we  derive  that,  for  any  x  in  IPJJ(A) 

III  G”x  III0.N  =  III  X  -  X  III0.N  <  III  X  -  N(N-l)<p(r(x))  1^  +  III  X  -  N(N- 1  )tp(r(x))  l^  ; 

hence 

III  (5^X  III0J4  <  (o(*Va‘)[  b(x  -  N(N-I)«|i(r(x),  X  -  N(N-l)(p(r(x))’'^ 

+  b(x-N(N-1)i|i(r(x),X-N(N-1)<p(r(x))’'^  1  . 
With  the  previous  notations  (see  (V.  1 4)  and  (V.  1 8))  and  recalling  that  the  normalization  factor 
between  the  forms  b  and  B  is  N(N-l)  (see  (111.14)),  It  is  simple  to  derive  that 
X  =  N(N-l)ip,^(r(x)).  while  X  (defined  in  (V.l  2))  satisfies  x  =  N(N-1  )ip,^/2(''(x));  we  derive 
III  (5^X  lllo,N  <  (oc*Vo(')(N(N-l))[  b(<jiN(r(x))  -  <p(r(x),  "Pn^Kx))  -  «p(r(x))’'^ 

+  b(  «h/2^^(x))  -  <p(r(x).  'Pn/2(''(x))  -  'P(r(x))’'^  ]  . 
applying  next  (V.20)  for  n  =  N  andn  =  N/2,  then  (V.15),  we  obtain 
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Lemma  V.2  :  There  exists  a  constant  c  independent  of  both  N  and  K  such  that 
(V.23)  V  X  €  IPj(A)  .  Ill  e"x  III0.N  <  clll  X  ill2.N  • 

Plugging  (V.23)  in  (V.2 1 )  and  using  (V.22)  give  now  the  following 

Theorem  V.  1  :  The  multigrid  algorithm  based  on  the  Jacobi  preconditioner  converges,  and  at  any 
V-cycle  with  m  smoothings  as  detailed  in  section  II  .  the  rate  of  convergence  is  given  by 
(V.24)  e;^(x)  dx  <  c /(m,1)  eQ^(x)  dx  . 

Remark  V.t  :  The  decay  rate  of  multigrid  algorithm  is  independent  of  N,  and  it  is  interesting  to 
note  that  its  asymptotic  behavior  is  0(  1  /m). 

Remark  V.2  :  The  analysis  of  the  case  of  non  constant  coefficients  requires  some  regularity  (i.e., 
)  of  the  coefficient  (  this  is  required  for  lemma  V.  1 ).  We  do  not  know  to  what  extend  the  lack 
of  regularity  of  oc  deteriorates  the  actual  convergence.  Note  however  that  the  local  regularity  is 
just  required,  i.e. ,  oc  can  be  non  smooth  through  the  Interfaces.  The  multigrid  procedure  proposed 
in  [2WH 1 .2  J,  and  that  is  based  on  another  approach  of  smoothing,  is  numerically  proved  to  be  very 
robust  with  respect  to  the  irregularity  of  oc  [*]. 

A.  APPENDIX 

The  main  purpose  of  this  appendix  is  to  provide  the  proof  of  lemma  V.  1 .  Let  us  first  recall 
the  following  result  of  the  approximation  theory  (see  [DJ) 

Lemma  A.  1  :  Let  n„  denote  the  orthogonal  projection  operator  from  L*(A)  onto  IP„(A).  The 
following  approximation  results  hold 

(A.  I )  V  u  €  H"’( A)  ,  B  u  -  n„,u  IIq  <  c(m)  M"*"  [  u‘’"’^(x)  ( 1  -x^)"*  dx]’^^  . 

Proof  :  Let  us  recall  that  the  Legendre  polynomials  constitute  a  total  system  of  orthogonal 

functions  of  L^(A)  that  verifies 

(A.2)  J,,L,(x)L„(x)(K  =  6„2/(2n.l)  , 


•  T.A.ZANG-  personnel  communication. 


i 
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Let  us  use  this  basis  to  span  the  function  u;  we  arrive  at 

u  L„  , 

so  that  the  L^(  A) -projection  of  u  onto  1Pm(A)  is  equal  to 

M 

“^n  :0  ' 

and  the  error 

(A  3)  u-nMU=2:„"M.,a„L„  • 

Let  us  recall  now  that  the  Legendre  polynomials  satisfy  the  following  relation 
(d-x^)  L„)'  =  -n(n+l)L„  . 

From  (A. 2)  we  conclude  that 

(A  "^)  (A'-„(x)L;^x)(1-x^)dx  =  8„^2n(n+l)/(2n+1)  . 

From  (A. 2)  and  (A. 3)  we  get 

llu- n^ullo  =  Lr--M.t  2a^/(2n+l)  , 
while,  from  (A.4) ,  we  derive  that 

u'^(x)  ( 1  -x^)  dx  =  2n(n+ 1 )  aj  /(2n+ 1 )  , 

and  (A.  1 )  is  then  just  a  simple  consequence  of  these  two  equalities  in  the  case  m  =  1 .  The  general 
case  is  handled  by  recursion. 

As  a  consequence,  we  derive  that 

Corollary  A.  1  :  Let  n,],  (tenote  the  orthogonal  projection  operator  from  H^CA)  onto  IP^CA). 
The  following  approximation  results  hold 

(A. 5)  Vu€H'^(A)nH’(A)  ,  ||  u  -  n^u  ||,  <  c(m)  M’"'  [  u‘"'*’’^(x)  ( 1 dx]’'^  . 

Proof ;  It  has  already  been  noted  that  for  any  u  in  Hq(A) 

[niu](x)  =  f*,  nM.,(u')(t)dt  , 
so  that  (A.5)  is  an  easy  consequence  of  lemma  A.  1 . 

Proof  of  lemma  V.1  :  It  is  an  easy  matter  to  note  that  the  solution  (ftp)  of  problem  (V.  1 7) 
satisfies 

Vx€A,  -  [(a[(p(p)],)J(x)  =  p(x)/(l-x^)  . 

From  (V.2)  and  (V.  16),  we  derive 

(A.6)  [<|)(p)lj(x)  dx  +  dx  <  c  p^(x)/(  1  -x^)  dx  <  c  B(p,e)  . 

From  corollary  A.  1 ,  we  know  that  there  exists  an  element  of  IPn/2(^^  • 
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I  liCe)  -  t'n  li  <  C  n’’  ['P(q)1J,(x)  ( 1  -x^)  dx  , 

and 

I  tnit  <1 'P(p)  It  . 
so  that 

(A.7)  I  i<)(p)  -  4.„  I,  ^  c  n’’  b(0,0)  . 

From  the  ellipticity  of  the  form  a" ,  we  derive  that  there  exists  a  constant  g,  such  that 

&  I  'Pn(e)  -  'Vn  iN  a'’('0n^p)  "  '*<„  .  'P„(p)  "  »J»„)  ; 

using  (V.  1 7)  and  (V.  1 8)  we  derive  that 

&  I  4>n(p)  -  tn  iN  "  '^r)  > 

or  again 

(A. 8)  &  I  (p„(0)  -  ‘i'n  I?  <  a('P(p)  -  'J'n  •  "PnCp)  “  »!'„)  +  SlC'I'n  *  "Pn^P^  “  '^r)  ~  ®"^'Pn  •  "Pn^P^  ~  'Pn^  • 

Let  US  note  now  that  from  the  exactness  of  the  Gauss  Lobatto  quadrature  formula  we  get  that 
(A, 9)  V  a„  €  (P„/2(A).  (<xy„  ,(4i„(p)  -  4-„)’)  =  (cx„<|/;  .(<P„(p)  -  . 

so  that 

a(<|(n  ,  'Pn^P)-<t'„)  -  a'’(4'n  .  'Pn(p)-<l'n)  =  ((a-0(„)i|»j,  ,(<?„( p)-4'n))  “  ((0C-a„)«|i^  ,(<P„(p)-'(>„))„ 
and  using  now  (A.8)  yields 

(A.  10)  I  <tln(p)  - 'I'nl,  <  c(  |(p(p)  -  +  l«  -  1  )  . 

It  is  standard  to  note  that  there  exists  an  element  a„  such  that 

n 


hence,  from  (A.  10)  and  (A.7),  we  derive  that 

I  -pCp)  -  -PnCp)  I,  <  c  n"’  (b(p.p)’'^  +  (Lj^o  I  '  'P^P^  li  ^  • 

and  (V.  1 9)  is  an  easy  consequence  of  (A.6). 

Let  us  turn  now  to  the  proof  of  (V.20).  We  shall  use  a  standard  duality  technique  and  first 
define  an  element  x  in  A)  as  follows 
(All)  Vv€Ho(A),  a(x,v)  =  b(((i(p)  -  t()„(p),  v)  . 

It  is  an  easy  matter  to  derive 

b(«p(p)  -  'p„(p).'p(p)  -  <Pn(p))  =  a(X.  <P(P)  -  <P„(P)) 

=  a(x  -  n’,2  X  .  'P(P)  -  <P„(P))  +  ^ 

=  a(x  -  Tt„’,2  X  .  <?(?)  -  «n^P))  +  (a"  -  aXftnVa  X  -  4>n(p))  • 

Let  us  examine  the  last  term  on  the  right-hand  side  of  this  last  equality;  using  (A.9)  one  more  time 
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we  write 

(a"  -  a)(n„’/2  X  .  =  ((«  -  ^n/2 

+  ((a-oc„)(  n„’^2  * 

so  that 

b(<p(p)  -  •Pn(Q).'p(e)  -  'Pn(e))  =  a(x  -  TI^2  X.1>(e)  -  'P„(e))  +  ««  -  an^(’^n/2  x)'.  'Pn(Q)')n 

+  ((a-oc„)(  n„’/2  x)‘.<P„(e)')  . 

hence 

b(<p(Q)  -  <p„(q).'p(p)  -  ipn^e))  < c|  X  -  n’/2^x I,  I (p(e)  -  ip„(e)  |, 

^  C  n-^  I  lr( A)  I  ’^n/2  X  li  I  «P„(e)  li  . 

<c|x-n^/2Xl,  l<p(p)-'Pn(p)li 

+  cn-2  2:J^ol«‘^’lL"•(A)l’^^2Xl,  l<P„(e)li  • 
Using  now  corollary  A.  1  and  (A.  1 1 ) ,  we  derive  (after  bounding  |  Il“*(  A)  ®  ^oi^stant) 

[b((p(Q)  -  (p„(p).ip(Q)  -  <P„(e))l^'^<c(n"’  |ip(p)  -ip„(p)|,  +  n'^|ip„(Q)|,  )  ; 
thanks  to  the  stability  of  (p^Cp)  with  respect  to  p,  we  derive  (V.20)  from  (Y.  1 9). 
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